clear all;clc;close all;
tic
cs=CubeSolver('nMax',11,'lambda',1.4,'nGrid',102,'lengthD',1.2,...
              'magD',0.2,'diffuseD',1.1,'fS',0.1,'fF',0.9,...
              'OD',2.2,'power',10.2,'P0',2.3,'tauD',1.2); % creat an object
toc

cs.dispEigvalue()
cs.disp1DCheck('cutOff',2)
cs.dispMode([1,2,1])
res=cs.eigenMode([1,2,2]);

w=0.2;
beamshape=@(x,y,z) exp(-((x-0.5).*(x-0.5) + (y-0.5).*(y-0.5))/2/w^2).*boolean(x>0.5);
res0=lambertw3D(cs, beamshape);

res1=cs.overlapMatrix(beamshape,'cutOff',2);
res2=cs.kappaMatrix('cutOff',2);
res3=cs.projectMatrix(beamshape,'cutOff',2);
cs.Solve( -res1 + res2, res3);
coeff=1.2;
res5=cs.lightPropagate(beamshape, coeff);

% test orthogonality
res6=cs.overlapMatrix(@(x,y,z) 1,'cutOff',2); % slow when cutOff is large

